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ABSTRACT 

We are concerned with the problem of buoyancy driven flow in a vertical, 
rectangular cavity whose vertical sides are at different temperatures and 
whose horizontal sides are insulated. An application of the dynamic A.D.I. 
method to obtain numerical solutions to this problem is described. For large 
non-dimensional temperature differences characterized by the Rayleigh number 
the flow patterns develop strong boundary layers. These boundary layers are 
resolved by applying the D. A.D.I. method to the discretization of this problem 
on a non-uniform grid. 
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Introduction 


The problem to be considered was proposed by Jones [4] as a test bed for 
numerical methods for solving a variety of practical thermal problems. It is 
known as the "double glazing" or "window cavity" problem and has many 
applications, particularly in the field of thermal insulation. The most well- 
known application is double glazing where a stagnant layer of air acts as an 
insulant between a warm room and a cold outside. The formulation of the 
problem is in terms of stream function, vorticity and temperature. The 
vorticity is eliminated to obtain a coupled pair of partial differential 
equations, one of which is fourth order. A dynamic A.D.I. (D.A.D.I.) method 
for solving these equations numerically is described. Grid stretching 
techniques originally due to Kalnay de Rivas [8] are used to resolve the 
boundary layers which develop for large values of the Rayleigh number. 


2. Formulation of the Problem 

We are concerned with the problem of fluid flow in an upright, 
rectangular cavity described in non-dimensional terms by 0 < x < 1, 

0 < z < 1, with z vertically upwards. The cavity has different constant 
temperatures on the vertical walls, T^ on the hot wall and T 2 on the cooler 
wall, and has insulated horizontal walls. We shall consider the two- 
dimensional flow of a Boussinesq fluid of Prandtl number 0.71 in which the 
flow takes place perpendicular to the walls. The Boussinesq approximation 
(see Mallinson and De Vahl Davis [9]) assumes that the physical properties and 
the density are constant except in the buoyancy term in the equations where 
the density is taken into account. This approximation is quite realistic and 
can give rise to predictions that are in good agreement with experiment (Jones 
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[5]) for small terraperature differences. The governing equations are 
considerably simplified by use of this approximation. A full discussion and 
detailed description of this problem is given by Jones [6]. 



Figure 1 . 


A non-dimensional temperature, T, is defined by 



where T* is the temperature. The equations representing the conservation of 
mass, momentum and energy may be written as 


V.v = 0, 

(Vxv) X v = - Vp - Ra Pr Tk - PrV 2 v, 


( 1 ) 


( 2 ) 
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V.(vT) = V 2 T, 


(3) 


where v = (u,0,w),k = (0,0,1) and p is the perturbation pressure. The 
Rayleigh number is given by 


Ra = ge(Tj-T 2 )/icv 


and the Prandtl number by 

Pr = v/k, 

where g is the acceleration due to gravity, 3 the coefficient of volumetric 
expansion, v the coefficient of kinematic viscosity and tc the coefficient 
of thermal conductivity. 

Mallinson and De Vahl Davis [10] show that the governing equations may be 
recast in the form 


? 3T Be 

Pr V £ + Pr Ra = u t— + w — , 
Bx 3x 9z 


= -C, 


S (uT) + sf (uT) ' ,2t - 


(4) 

(5) 

( 6 ) 


u 


3j/ 

3z 



(7) 


where ? is known as the vorticity and the stream function. Eqs. (4) to 
(7) represent what is known as the vorticity - stream function formulation of 
the problem. 
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The boundary conditions are 

* = H = °» on r « 

where T is the boundary of ft = {(x,z): 0< x < 1,0< z< 1}, 

T = 1, on x = 0, 

T = 0, on x = 1, 

3T 

“ 0, on z = 0 and z - 1. 

Eliminating the velocities u,w and the vorticity £ we obtain the following 

system of equations: 




Pr 3x 3z ^ V) 


= 0 , 


( 8 ) 


V 2 T - («T) - £ («T) - 0. 


(9) 


An advantage of eliminating the vorticity is that the need for a vorticity 
boundary condition is avoided. The system of equations (8) and (9) represents 
a fourth order equation for ^ and a second order equation for T. The major 
interest is in heat transfer so there are two important quantities namely the 
temperature fields and the overall heat transfer defined by a Nusselt number 
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In order to solve this problem using an A.D.I. method we must first convert 
Eqs. (8) and (9) to the parabolic equations 


I 

X 3t 


vS - 



+ _L |1 _1 r v 2^ _ _i li _J_ r v 2^ „ 0 

Pr 3z 3x 1 Pr 3x 3z 1 VJ 


( 10 ) 


3T 

3t 



(uT) 



(ID 


whose steady state solution, if one exists, solves Eqs. (8) and (9). The 
parameter X controls the interaction between the equations. It means that 
for X * 1 we are effectively using different time scales for the two 
equations . 


3. Solution of a Nonlinear Equation 

In a numerical process to find the solution of Eqs. (8) and (9) we use 
Eq. (8) to solve for ip and Eq. (9) to solve for T. Eq. (8) is nonlinear 
with respect to ip and can be written in the form 

L (ip) = f (T) , (12) 

where L is a nonlinear operator. 

* 

A Newton-type method is used to solve Eq. (12). Suppose that ip is 

some approximation to the solution of Eq. (12). We replace L by its 

* 

linearization about ip and then attempt to partially solve the linearized 


problem 


L + lO*) - f (T) = 0 


(13) 



6 


where L ) is the Frechet derivative of L at ^ . We use the D.A.D.I. 
method to solve Eq. (13), the linearization being updated after each D.A.D.I. 
step. This will be explained more fully later. 

The Frechet derivative of L at is given by 


L'rt)-* = V% +? i |t 


Pr L 3x^ 


- si 


Pr 1 3 z^ 


3x 


(14) 


4. Finite Difference Equations 

We cover ft with a square grid of mesh length h = 1/N where N is a 
positive integer. Let , and T be the values of and T at the 

i y J 1 > J 

point (x^,Zj) respectively where = ih and z^ = jh. We use standard 

second-order centered difference approximations. At this point it is 

2 — 

convenient to introduce the central difference operators 6 and 5 as 

x x 


follows: 


5 v . , = v . . , . - 2v . + v 4 t , 6 v . . = v . , - . - v . - 

x i, j i+l,j i»j i-l»j x i,j i+l,j i-l»; 


The central difference operators 6^ and 6^ are defined in a similar 
manner. Eqs. (8) and (9) are discretized using these approximations to give 


-4 (S 4 +26 2 6 2 +6 4 U , -nk~ (u. ,6 p. . + w. .6 p .) = 0,(15) 

4 ^ x z x z jr l,3 2h x i, j 2hPr i,j x p i,j i,j z^i ,j J 


and 


J - 2h V“ T >l,j - 2h S |“ T >i,J " °- 


(16) 
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where 


u 


ij 


“ 2h 5 6 z*±,y W ±,j 2h 6 x*i,J* P i,j .2 


The normal derivative boundary conditions are also discretized using 

central differences and these are used to remove imaginary exterior points 

from Eqs. (15) and (16) near the boundary. We have 2N(N-1) equations for 

the unknowns ,i = 1,...,N-1, j = 1,...,N-1, and T ,i = 1,...,N-1, 

i>j 

j = 0 , • • • ,N • 


Eq. (15) can be written in the form 






which is the discrete form of Eq. (12), where is a nonlinear discrete 

operator. 


5. Method of Solution 

Here we describe how the D.A.D.I. method of Doss and Miller [2] can be 
applied to this problem. Since Eqs. (10) and (11) are parabolic they may be 
advanced in time by a direct method and the complete solution procedure may be 
regarded as a single iterative scheme. Our interest is not in solving the 
parabolic equations (10) and (11) accurately for finite times but to reach the 
steady state solution as soon as possible. We therefore use the D.A.D.I. 
method which uses a strategy that attempts to keep the time step At within a 
region of fast convergence. An advantage of using an automatic step size 
changer is that it avoids the necessity of choosing a priori iteration 
parameters. The strategy of Doss and Miller [2] attempts to recognize 
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instabilites as they start to occur and to bypass them by decreasing At. 

We start with initial approximations and to T and ^ 

respectively at time t - 0. Suppose that we have reached the nth time step, 
where n = 2m and m is even, and that the current approximations of T and 
ij> at the grid points are and respectively. Starting from these 

approximations we describe a step of the D.A.D.I. method with time step 
At. 

We begin by putting 


for i , j “ 0,1, ...,N. 

The A.D.I. scheme due to Peaceman and Rachford [11] is used to advance 
Eq. (11) in time. The following systems are solved along lines in the x- 
direction: 




2v(n) 


i = 


(17) 


(r^-S 2 )^ + V 2 h6 x (u (tl) T< a+1) ) 1> . 


= - 1 /2^ z (w (n) T (n) ) i>j 


(18) 


i = 1,...,N-1, j = 


with 


= ( r. 1+6^ )T^ n ^ 
v 1 x' i,N ^ 1 z J i,N 


2w(n) 


(n+1) _ . (n+1) _ - 

T 0,j " T N,j " °» 


i = 1, 


j 0, 1 * • • • »N, 


( 19 ) 
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homogeneous Neumann boundary conditions along z = 0 and z = 1 and where 
r^ = h^/At. 

The following systems are solved along lines in the z- direction: 


(V-'MrT +1 /2h«> (n) T (n+2) ) lid 


- (r,I+6 2 )T <n t 1) -V 2 h? (u (n) T (n+1) l 4 

1 x' i,j * x" i , j 


( 20 ) 


j = 0, * . . ,N, i = 


with homogeneous Neumann boundary conditions along z = 0 and z = 1. 

We advance ^ in time using the A.D.I. scheme of Douglas and Rachford 
[3] to solve the fourth order linear equation (13). The extension of this 
scheme to solve the biharmonic equation is due to Conte and Dames [1]. 

We solve the following systems, described by j = along lines 

in the x- direction: 


- <V-<-2A 2 - 457 [* + 457 


- 457 + «) - 


(n) 


,(n)> 


i = 1,...,N-1, 


( 21 ) 


with homogeneous Dirichlet and Neumann boundary conditions, and where 
r 2 = h^/(XAt). 
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The following equations, described by i = 1,...,N-1, are solved along 
lines in the z- direction: 


{r,I+S 4 - 75 - (6 i|> (n })6 6 2 +-t|- [6 (6 2 +5 2 )^ (n J ]fi }$J n 'J' 2) 
1 2 z 4Pr v x i,j' z z 4Pr L x^ x z J i,j J z JT i,j 


t 5 '- ( 5 K+i?>z S z + 4fe + '2*^ U - (22) 


j = 


We define an A.D.I. step to be the process by which we obtain T 


(n+ 2 ) 


and ,j)( n+ 2) from T^ and <j>^ i.e. the solution of equations (17) to 
(22). Starting with the approximations T^ n+ ^ and <f,( n+ 2) we perform a 


second A.D.I. step using the same time step At to obtain T 


(n+4) 


and 


<))(n +4 ). The next part of the process is the bookkeeping stage of the 

D. A.D.I. step. Here we start with the approximations and <|>^ n ^ and 

( n+4 ^ ^ ( n+4 ^ 

perform an A.D.I. step with time step 2At to obtain T and <j> . 

We compute the test parameter, TP, which is given by 


TP = / [SUM/ASUM] , 


where 


and 


SUM = llT (n+4) -? (n+4) H 2 + ll<}> (n+4) 4 (n+4) ll 2 


ASUM = HT (n+4) -T (n) ll 2 + «*< n+4 >-* (n) | 2 . 


The test parameter is an estimate of the relative local truncation error. If 
we are interested in solving the parabolic equations ( 10 ) and ( 11 ) accurately 
then At will be small, and so will TP. Our main concern, however, is to 
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accelerate convergence and attempt to push TP into an interval where 
convergence for A.D.I. is rapid* The strategy for changing At is described 
in Doss and Miller [2]. Although their analysis of the step size strategy 
rests on rigid assumptions, they obtain good results in situations beyond the 
realm of the assumptions* If TP > 0.6 then we reject the present D.A.D.I. 
step and start the step again with At reduced by a factor of 1/16* If 
TP < 0.6 then we accept the present D.A.D.I. step and change At according 
to the strategy: if TP falls in the intervals [0,0.05], (0.05,0.1], 

(0.1, 0.3], (0.3, 0.4], (0.4, 0.6] change At by the factors 4,2,/3, V2 > V4 

respectively for the next D.A.D.I. step. 

If the present step is accepted then we use one step of Newton's method 
to update the approximation to by setting 

^(n+4) = # (n> _ ^ (n+4 ) ^ 

We are now in a position to start the next D.A.D.I. step. 


6. Numerical Results 

The initial approximations T^ 0 ^ and are chosen to be the values 

at the grid points of the functions T^ and ^ respectively, where 

Tj (x,y) = 1 - x, ^ T (x,y) = 0. 

The algorithm described in the previous section was run for different grid 
sizes and for Rayleigh numbers of 10^,10^,10^, and 10^. The algorithm is 


terminated when the maximum modulus of the difference between successive 
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iterates on even time steps is less than 10"*^ for h = 1/16,1/32 and 
5 x 1(T 5 for h = 1/64. 

We experimented with different values of the parameter X. For X = 1 
we found that for large values of the Rayleigh number the iterates oscillated 
due to the interaction between the equations. For suitable choices of this 
parameter the instabilities were controlled and the method converged. The 
value of X was chosen to be 0.05, 0.01, 0.002 and 0.002 respectively for 
Rayleigh numbers of 10^,10^,10^ and 10^. The choice of X was not 

critical in the sense that similar convergence behavior was observed for a 
wide range of values in the neighborhood of the chosen value. 

Suppose that for a particular value of the Rayleigh number the problem 
has been solved numerically on a grid of mesh size h. To find the solution 
on a grid of mesh size V 2 h we use the best available information to begin 
the new calculation i.e* we use as our initial approximation values 
interpolated from those obtained on the coarser grid. Cubic interpolation is 
used for values of and linear interpolation for values of T. 

The average value of the Nusselt number is calculated using the 
trapezoidal rule, where we use the following approximation to the normal 
derivative of T on the hot wall: 

3T,„ , .. {-3T(0,z)+4T(h,z)-T(2h,z)} 

3^ ( °’ Z) " 2h 

Table I contains the average Nusselt number on the hot wall and the 
maximum and minimum local Nusselt numbers on the hot wall, and their 
location. Table II contains the number of D.A.D.I. steps and the run time, in 
seconds, required to reach the convergence criterion. The results were 
obtained on the Oxford University ICL 2980 computer. Contours for the 
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temperature and stream function are shown in Figures 2 to 9 for different 
values of the Rayleigh number. 

The results we have obtained mainly show a good agreement and the 
convergence is approximately quadratic. Also it is clear that some form of 
non-uniform grid is needed to resolve the boundary layers accurately when 
Ra is large. In the next section we show that this technique has the effect 
of improving the rate of D.A.D.I. convergence. 


TABLE I 

Heat Transfer Results 






Average 

Maximum 

Minimum 

Ra 

Mesh 

Nusselt Number 

Nusselt 

Number 

Nusselt 

Number 





Nu 

^ u max 

@z = 

Nu min 

@z = 


17 

X 

17 

1.120 

1.515 

0.86 

0.697 

0.00 

10 3 

33 

X 

33 

1.118 

1.508 

0.91 

0.693 

0.00 


65 

X 

65 

1.118 

1.507 

0.91 

0.692 

0.00 


17 

X 

17 

2.357 

3.838 

0.81 

0.615 

0.00 

10 4 

33 

X 

33 

2.270 

3.628 

0.84 

0.592 

0.00 


65 

X 

65 

2.250 

3.554 

0.86 

0.587 

0.00 


17 x 

17 

5.146 

8.423 

0.88 

0.854 

0.00 

33 x 

33 

4.758 

8.508 

0.91 

0.755 

0.00 

65 x 

65 

4.573 

7.972 

0.92 

0.735 

0.00 


lCr 


33 x 33 
65 x 65 


10.062 19.085 0.91 0.926 

9.272 19.590 0.95 0.999 


0.03 

0.00 


x 
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3 

Figure 2. Temperature contours for Ra = 10 . 
Contour levels are 0.1(0. 1)0.9. 



Figure 3. Stream function contours for Ra = 10^. 

Contour levels are -0. 12 (-0. 1 2 ) — 1 .08. 
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Figure 6. Temperature contours for Ra = 10^. 

Contour levels are 0.1(0. 1)0.9. 



Figure 7. Stream function contours for Ra = 10^. 

Contour levels are -1.06(-1.06)-9.54. 
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7. Non-Uniform Grids 

For large values of the Rayleigh number the flows develop strong boundary 
layers. This effect can be seen in Figure 9, for example. When Ra = 10^ we 
have had to use an extremely fine mesh, 65 x 65, in order to resolve these 
boundary layers. This is rather wasteful since the grid points are also 
densely distributed away from these layers where they are not needed. Here we 
apply the technique of Kalnay de Rivas [8] to resolve the boundary layers 
using a non-uniform grid. 

Basically, the idea is to make a change of independent variable so that 
the domain is mapped into a new co-ordinate system where the variations of the 
solution are not so rapid. The grid intervals are varied by defining 

stretched co-ordinates C and n such that x = x(?) and z = z(n) where 
the grid intervals A? and An are constant and x and z are the old 

physical co-ordinates. The mapping is chosen so that the solution, when 

regarded as a function of the new variables, has no boundary layers. 

Kalnay de Rivas [8] and Jones and Thompson [7] show how to express 

derivatives in terras of the stretched co-ordianates. For example, we can 
express the first derivative in terms of ? in the following manner 


8 v _ 8v # d£ 
3x 8£ dx* 


(23) 


Eq. (23) can be discretized using central differences to give the 
following approximation 


9 v B V l+l.i~ V l-1.1 d£| 

8x 2A£ dx|x^’ 


where v^ j is the value of v at the grid point (x^,Zj) with 


( 24 ) 
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x^ = x(iAc) and z^ = z(jAn). 

The transformation can be differentiated using central differences to 
obtain the following approximation to the first derivative: 

9x 2AC ^i+r x i-i^ i x 1+ i- x i_ii 

Finite difference approximations for higher order derivatives are obtained in 
a similar way. 

Let x(C) and z(n) be two grid stretching functions with constant grid 

intervals A£ and An respectively. The region Q is covered with a 

variable grid defined by the above mappings. We define T. , and . to 

- L >J i> j 

be the values of T(x,z) and ij>(x,z) respectively at the grid point 
(x^,Zj). The finite difference equations are constructed using approximations 
like Eq. (25) for the derivatives. The system of finite difference equations 
are solved using the D.A.D.I. method described earlier. 

The problem is solved on the following non-uniform grids: 


(a) x(C) = sin^( V 2 nc) , 

z(n) = sin ( V 2 Hn) • 

(b) x(C) = 6? 5 - 15c 4 + IOC 3 , 

z(n) = 6n 5 - 15n 4 + 10n 3 . 


(26) 


(27) 


These grid stretching functions give a smaller spacing of the grid points near 
the boundary. 
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The problem was solved for Rayleigh numbers of 10^ and 10^ on the non- 

uniform grids defined by Eqs. (26) and (27) with A£ = An = 0.04. In Table 

III we compare the average value of the Nusselt number along the hot wall 
obtained on these non-uniform grids, a uniform 26 * 26 mesh and a uniform 

65 x 65 mesh* In Tables IV and V we compare the number of D.A.D.I. steps 

and run time respectively to reach the convergence criterion, which is the 
same as that used for the results in Table II. 

In Table III we see that we have obtained a good estimate of the Nusselt 
number by using a stretched grid even though we have only a 26 x 26 mesh. 
For small values of the Rayleigh number the use of a stretched grid has little 
effect. It is only when the Rayleigh number becomes large that we obtain an 
improvement by using non-uniform grids. The use of non-uniform grids to 
resolve the boundary layers has had the effect of improving the rate of 
D.A.D.I. convergence. Approximately twice as many steps of D.A.D.I. were 
required for convergence of the method using the uniform 26 x 26 mesh than 
for stretched 26 x 26 meshes. 

The use of non-uniform grids has allowed us to resolve the boundary 
layers using comparatively few mesh points. Reasonable accuracy has also been 
obtained on the stretched grids compared with an extremely fine mesh. 

The iterative solution of finite difference equations constructed on a 
non-uniform grid usually presents great difficulties. This is due to the 
problem of finding suitable parameters for the acceleration of convergence of 
any selected iterative method. Hence, an advantage of the D.A.D.I. method 


over standard iterative methods for solving problems of this type is that we 
do not require an a priori choice of parameters to accelerate convergence. 
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TABLE III 

Average value of the Nusselt number 



Uniform Mesh 

Stretched Mesh 

Stretched Mesh 

Uniform Mesh 

Ra 


(a) 

(b) 



26 x 26 

26 x 26 

26 x 26 

65 x 65 

10 5 

4.889 

4.596 

4.595 

4.573 

10 6 

10.163 

9.123 

9.066 

9.272 




TABLE IV 

Number of D.A.D.I. 

steps 



Uniform Mesh 

Stretched Mesh 

Stretched Mesh 

Uniform Mesh 

Ra 


(a) 

(b) 



26 x 26 

26 x 26 

26 x 26 

65 x 65 

10 5 

151 

85 

84 

80 

10 6 

676 

334 

369 

426 




TABLE V 
Computational 

time 



Uniform Mesh 

Stretched Mesh 

Stretched Mesh 

Uniform Mesh 

Ra 


(a) 

(b) 



26 x 26 

26 x 26 

26 x 26 

65 x 65 

10 5 

68 

85 

84 

238 

10 6 

299 

332 

364 

1280 
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